Simulations of strongly phase-separated liquid-gas systems 
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Lattice Boltzmann simulations of liquid-gas systems are believed to be restricted to modest density 
ratios of less than 10 In this article we show that reducing the speed of sound and, just as 
importantly, the interfacial contributions to the pressure allows lattice Boltzmann simulations to 
achieve high density ratios of 1000 or more. We also present explicit expressions for the limits of 
the parameter region in which the method gives accurate results. There are two separate limiting 
phenomena. The first is the stability of the bulk liquid phase. This consideration is specific to 
lattice Boltzmann methods. The second is a general argument for the interface discretization that 
applies to any diffuse interface method. 



Simulations of liquid gas systems with lattice Boltz- 
mann have been restricted to small density ratios in the 
past. Those restrictions have lead to the development of 
hybrid lattice Boltzmann schemes to be able to simulate 
systems with high density ratios of about 100-1000 by 
Inamuro et a/.|l(. I n this article we explain how large 
density ratios can also be achieved with standard lattice 
Boltzmann methods. Furthermore we derive the condi- 
tions which limit the ability of the method to obtain sta- 
ble, accurate, and unique results for the phase diagram. 
We present a new general argument for the minimum in- 
terface width required to accurately simulate a system at 
a given reduced temperature. This important argument 
is not restricted to lattice Boltzmann methods but only 
relies on the relation of the discretized interface to the 
expression for the pressure. It therefore applies to all 
diffuse interface methods. 

The lattice Boltzmann method can be viewed as a dis- 
cretization of the Boltzmann equation. The hydrody- 
namic limit of the Boltzmann equation gives the conti- 
nuity and Navier-Stokes equations and the discretization 
of the lattice version is chosen such that it preserves this 
limit. The basic variables of the lattice Boltzmann equa- 
tion are a set of densities /i(x, t) associated with a ve- 
locity set Vj. The evolution equation for the fi is then 
given by @ 

/i(x+Vi,t + l) = / < (x,t) + F < + i(/P(x,t)+^ i -/ i (x,t)). 

T 

(1) 

The are the equilibrium distribution corresponding 
to the ideal gas. Non-ideal contributions are included 
through the bulk forcing term Fi or pressure term Ai, 
following 0- The fluid density is denned as p = Ei/>) 
and the momentum is pu = fiVi (although the to- 
tal momentum contains additional contributions from the 
force) . The moments of the equilibrium distribution are 

Ei f° = P, E* fi v ia = pu a , 
Ei fiViciVip = pu a U/3 + p95 al3 , 
Ei fi v iaVipVi 7 = p9(u a 8f3 7 + Ufs5 al + u 7 S a/3 ) 

+pu a U l3 U 1 + Qccfr- (2) 



Here Q is a correction term that should be zero. Most 
velocity sets for lattice Boltzmann are limited to Vi X S 
{ — 1,0, 1} so that vf x = Vi X . This restricts the third mo- 
ment in J2J) to 9 = 1/3 and Q a /3-y = —pu a upu-y. 

The non-ideal contributions from the Ai need to con- 
serve mass and momentum and the moments are given 
by 

E^ = 0, E<^(vi-u) = 0; 

J2i Ai{v ia - up)(vip -up) = A a p, (3) 

Ej Ai{v ia - U a ){Vip - Uf3)(v i7 - Uj) = AapUy + A ai up 

+ A/3yU a . 

The forcing term Fi has the moments 

Ei^=0, Ei^(vi-u)=F, 

EF i (v i -u)(v i -u) = *. (4) 

A standard expansion of (yi gives the continuity equa- 
tion 

d t p + V(pu)=0 (5) 

and the Navier Stokes equation 

dt(pu) + V(puu) = -V(p0 + A) + F + Vct + VR (6) 

where u = u + ^F 0- Here the Newtonian stress tensor 
is given by 

a = V p{Vu + (Vu) T ) (7) 

and unphysical terms have been collected in the remain- 
der tensor 

R = t* - 3v[uV.A + (uV.A) T + u.VAl + VQ] + O{0 2 ). 

(8) 

The kinematic viscosity is given by v = (t - \)9. Note 
that most of the unphysical terms in JHJ violate Galilean 
invariance 0,0. 

We see from (jHJl that for A = and F = the lattice 
Boltzmann method enforces an ideal gas equation of state 
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FIG. 1: Sound speed squared c s for the liquid and gas phases 
as a function of the quench depth for po — 1. 



with p(p, 6) = p9 = p/3. To simulate a fluid with a non- 
ideal equation of state P(p, 9) — p9 + P md (p, 9) we can 
now choose 



F = — V.P 

* = (T-\)pFF+±(X7V) D p, 
A = 0, 



(9) 



which we will refer to as the forcing method |2| • The care- 
ful reader may have recognized that the W term does not 
appear in the Navier-Stokes equation ©, but a higher 
order analysis shows that these terms are necessary to 
recover the correct equilibrium behavior 2]- An alterna- 
tive choice for the moments is 



F 

A : 



0, W 

pnid | 



= o, 

v{uS7 p 



(uWpf + u.W P 91), (10) 

which we will call the pressure method 0, Q . For either 
approach we recover the Navier Stokes equation for a 
non-ideal gas 



d t (pu) + V(puu) = -VP + V(T. 



(11) 



In equilibrium both approaches lead to a constant pres- 
sure P and therefore to the same density profiles Q. 

Most previous lattice Boltzmann simulations ap- 
proached the simulation of non-ideal systems by using 
the ideal gas equation of state p = p9 = p/3, as a starting 
point. Interactions are then included to allow the simula- 
tion of non-ideal systems. The speed of sound c s = U d p p 
will then recover the ideal gas value of 1/3 in the dilute 
limit. For a van der Waals gas with a critical density of 
1 and a temperature of 9 — 1/3 and an interfacial free 
energy of J ft/2 (Vp) 2 the pressure tensor is given by 



P 



Po 



3-p 



IP 2 



-K [ pV z p + -Vp.Vp ) 1 + nVpVp 



(12) 
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FIG. 2: The recovery of the phase diagram for different values 
of the interface width. The van der Waals phase diagram is re- 
covered to very good approximation. For very deep quenches 
corresponding to large density ratios wide interfaces are re- 
quired to recover the very low gas densities. The value of po 
does not affect the form of the interface or the value of the 
gas and liquid densities and values from 1 to 10 -7 were used 
for increasing quench depth. 



Previous approaches matched the ideal gas equation of 
state in the dilute limit, leading to po = 1. For the 
van der Waals gas the speed of sound increases rapidly 
for high densities. A problem arises when the speed of 
sound becomes larger then the lattice velocity \vi\, be- 
cause information can not be passed on at speeds larger 
than the lattice velocity. When the speed of sound is in- 
creased above 1 the simulation becomes unstable. This 
clearly limits the range of critical temperatures for which 
we can obtain stable solutions in lattice Boltzmann, as 
shown in Figure ^ 

The stability analysis is slightly complicated by the 
fact that we have additional gradient terms in the pres- 
sure tensor. These terms further decrease the stability, 
as shown in a previous analysis of the pressure method 
by C. Pooley for one, two, and three dimensional lattice 
Boltzmann methods 5] . In the notation of this letter the 
linear stability condition is 



(13) 



for a homogeneous system with density p. This suggests 
that, at least as far as the stability of the bulk phase is 
concerned, the most stable solutions should be found for 
k = 0. 

To lower the speed of sound in the liquid phase we now 
reduce the value of po in ltl"2"|l . This decreases the speed 
of sound in the liquid by a factor pq. This also increases 
the range of stability for n in (|13|l . We now expect that 
lowering the speed of sound by a sufficient factor will 
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(a) Pressure Method 




(b) Forcing Method 



FIG. 3: Existence of accurate solutions for the (a) pressure 
and (b) corrected forcing method for different values of po 
and w. Symbols indicate parameter combinations that lead 
to stable, accurate, and unique solutions. Solid lines are the 
bulk stability limits for the pressure method given by eqn. 
1131 . The dashed line is the line for an accurate interface 
representation given by eqn. 1171 . 



reduce the speed of sound sufficiently to simulate systems 
with arbitrarily low temperature ratios 



To test this idea we performed simulations with near 
equilibrium profiles using a one dimensional three veloc- 



ity Vi = { — 1, 0, 1} model by denning 
Pi - Pa , 
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where pi and p g are the equilibrium gas and liquid densi- 
ties and N x is the number of lattice sites. The interface 
width is given by 



w(k,0/0 c ) = 



2 k 



%/0 - l 



(15) 



This profile is not the exact analytical solution to the 
differential equation VP = 0, but it is very close to it. 
By initializing the simulation with this profile we can test 
the linear stability of the method around an equilibrium 
profile to good accuracy. The shape of a stable interfacial 
profile is independent of po for both the pressure and the 
forcing method. 

In Figure [5] we see that by lowering po the method is 
now able so simulate very small values of the reduced 
temperature 0/0 c for interface width w > 1, but that 
significantly larger width are required to recover an ac- 
curate phase diagram for deep quenches. For values of 
0/0 c between 0.9 and 1 we also find non unique solutions 
for small values of k which is discussed in more detail in 
a previous paper |2(. 

To understand when the method fails to obtain ac- 
curate results note that, in equilibrium, equation 112|) 
requires that 

_ p/(3 -p)- 9/8 P H 



Pi, 



pS7 2 p- iVp.Vp 



(16) 



where pb is the bulk pressure corresponding to the den- 
sity p. For small values of k the interface becomes sharp 
in the continuous limit so that the derivatives become ar- 
bitrarily large. But in the numerical implementation the 
derivatives are discrete derivatives. The discrete values 
are limited by the lattice spacing. In the one dimen- 
sional case we choose Vp(x) = 0.5(p(x + 1) — p(x — 1)) 
and V 2 p(x) = p(x + 1) — 2p(x) + p(x — 1). For higher 
dimensional stencils with the same stability limits for the 
bulk phase see C. Pooley's thesis 0- 

The methods always lead to a constant pressure, even 
across an interface^. We can now perform a simple esti- 
mate of the minimum value K m that allows this pressure 
to be the equilibrium pressure. For any point with den- 
sity p s we can consider two neighboring points, one with 
a smaller density p_ and one with a larger density p+. 
We can now find a lower limit for the smallest value K m 
by varying the values of p+ and p_ 



max 

Pg<Ps<Pl 



mm 

Pl<P-<p s 
Ps<P+<Pl 



p/(3 -p)- 9/8 p 2 



Pb 



t(j>- - 2 Ps + P+) - \{P-\ 



(17) 
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where pi is the liquid density and p g is the gas density. 
We performed a scan of the parameter space w and 0/9 c 
initializing the simulation with a near equilibrium profile 
for different values of po- We will accept simulations that 
are stable, accurate and unique. We choose as the crite- 
rion of accuracy that log 10 (/5 m m) — l°gio(/°g) < 0.1. As 
can be seen in Figure |2J the results are not very sensitive 
to the exact value of the cutoff. For values of the inter- 
face width w < 1.5 we also test the uniqueness of the 
simulation by using initial profiles with bulk densities 
corresponding to the pressure at the spinodal points^- 
Our criterion for uniqueness is then that all simulations 
lead to the same minimum density to within Ap < 0.01. 

The comparing (|17|l . shown as a dashed line in Fig- 
ure 01 and the numerical results for stable, accurate and 
unique solutions shows excellent agreement. The bulk 
stability of eqn. 113fl gives the second limit for the ac- 
ceptable parameter range for the pressure method. The 
forcing method leads to a slightly larger range of bulk sta- 
bility. The underlying reason is that calculating deriva- 
tives in Q leads to an additional information exchange 
allowing for speeds of sound slightly larger than 1. But 
the dependence of the stability on po and w is very simi- 
lar to the one for the pressure method. Note that previ- 
ous lattice Boltzmann simulations correspond to po = 
which corresponds to a small acceptable parameter range. 

The interface constraint (|17|) is remarkably successful 
at predicting the acceptable simulation parameters. It 
predicts how thin is too thin for an interface. It thereby 
detects when non-unique solutions occur and when solu- 
tions for deep quenches fail to deliver accurate results. 
The criterion presented so far is entirely numerical but 
because of its importance we want to examine two lim- 
iting cases for shallow and deep quenches for which we 
can obtain analytical results. 

We first examine for which values of p s the minimum 
value of k is reached in eqn. I|17[l . The dashed line in the 
inset of Figure 0] shows how this density p cr a varies as a 
function of temperature. 

Near to the critical temperature, the orange line in the 
inset in Figure 01 lies close to the high density spinodal 
curve. This is because P — pb has its highest magnitude 
here, therefore helping to maximize n m within this re- 
gion. An analytical estimate for be obtained by 
expanding the pressure around the critical density, giving 
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a Analytic Solution Low 6/9. 

□ Analytic Solution High 6/9 



10 



FIG. 4: The two limiting cases for which we can obtain an 
analytical approximation to the w{K m ) relation. The inset 
shows the value of the critical density p cr it for which the most 
severe limitation for K m CZjl occurs. Note that there is a 
discontinuity. 



into l|15[) gives a minimum interface width of 



(19) 



As the temperature is decreased in the inset of Figure 0] 
the critical density p cr n makes a discontinuous jump to 
a regime in which it lies close to the gas density, p g . The 
minimum interface width can, in this case, be analytically 
obtained by expanding densities around p g . We define 
p s = p g + dp and p + = p_ + Ap. Since P — pb is a 
positive quantity, a suitable choice for p_ is p_ = p a . 
Substituting these expressions into equation 1|17[) gives 



(20) 



(Pg + 8p)Ap- Ap 2 / 



Minimizing this with respect to Ap leads to p s — Ap/4. 
Re-substituting this result back into equation 1)20(1 . and 
maximizing with respect to 5, we finally obtain p cr it = 
2p g . Using this we can calculate the minimum interface 
width, 



P 



Pb 



0) (Ps -l) + -(p s - if ■ (18) 



Within this regime, P—pb is large and negative, therefore 
p_ and p+ must be chosen to make the denominator in 
()17J) as negative as possible. A suitable choice is p_ = p g 
and p+ = p s . We assume that the critical value of p s lies 
on the spinodal curve p sp i n — 1 + 2\/0 c — 6. This allows 
us to obtain K m (p cr it), and substituting this expression 



1 



V^Pg Wo ~ 0) 



(21) 



as shown by the triangles in Figure^] This closely follows 
the numerical result at low temperatures. 

We have therefore shown how to simulate deep 
quenches with lattice Boltzmann and we were able to 
predict which simulation parameters will lead to accu- 
rate results. 
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